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Inspiraling supermassive black hole binary systems with high orbital eccentricity are important 
sources for space-based gravitational wave (GW) observatories like the Laser Interferometer Space 
Antenna (LISA). Eccentricity adds orbital harmonics to the Fourier transform of the GW signal, 
and relativistic pericenter precession leads to a three-way splitting of each harmonic peak. We study 
the parameter estimation accuracy for such waveforms with different initial eccentricity using the 
Fisher matrix method and a Monte Carlo sampling of the initial binary orientation. The eccen- 
tricity improves the parameter estimation by breaking degeneracies between different parameters. 
In particular, we find that the source localization precision improves significantly for higher-mass 
binaries due to eccentricity. The typical sky position errors are ~ 1 deg for a nonspinning, 10 7 M© 
equal-mass binary at redshift z = 1, if the initial eccentricity 1 yr before merger is eo ~ 0.6. Peri- 
center precession does not affect the source localization accuracy significantly, but it does further 
improve the mass and eccentricity estimation accuracy systematically by a factor of 3-10 for masses 
between 10 6 M© and 10 7 Af© for e ~ 0.3. 

PACS numbers: 04.30.Db, 04.80.Nn, 97.60.Lf 
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I. INTRODUCTION 

The inspiral and merger of compact binary systems of 
black holes are important sources of gravitational waves 
(GWs) for the proposed space-based GW missions such 
as the Laser Interferometer Space Antenna (LISA) [l[ 
or the European New Gravitational Wave Observatory 
(NGO/eLISA) gj. The detectable frequency band for 
these instruments will be around 1CU 4 to 10 _1 Hz Q 
which corresponds to the inspiral of two (10 4 — 10 7 )M Q 
black holes. As the sources detected by LISA/NGO will 
be loud with a large signal-to-noise ratio in general, an 
ideal method for parameter extraction is matched filter- 
ing 0. 

An effective matched filtering requires an accurate 
model of the emitted GWs. In this technique the de- 
tected signal output is cross-correlated with theoretical 
waveform templates. In particular, matched filtering is 
sensitive to the phase information of the waveform, and a 
high correlation between the signal and template allows 
one to make predictions on the source parameters 0, Q . 

Many previous studies in the literature adopted wave- 
forms generated by binaries in circular orbits (see Refs. 
0-Qjl for LISA parameter estimation). This is due to the 
expectation that the orbit of the binary will circularize 
due to the emission of GWs [H |2l|. 

Nevertheless, there are a number of reasons to expect 
that at least some LISA sources may be eccentric. If the 
binary is embedded in a gaseous disk, it can remain ec- 
centric until the final year of the inspiral [22H25I ] . The 
interaction of the supermassive black hole (SMBH) bi- 
nary with a population of stars also increases its eccen- 
tricity [2a - |28| . The eccentricity can also be excited by 
the Kozai mechanism and relativistic orbital resonances 



in hierarchial triples [29M33l | or by a triaxial potential 
[34l l35ll . and may be typical for extreme mass ratio inspi- 
rals |36l. l37j . Further, black hole binaries in dense galac- 
tic nuclei formed by GW emission during close encoun- 
ters remain very eccentric until merger (38l . |39| . Popula- 
tion synthesis and binary evolutionary models show that 
a fraction of stellar compact object binaries may also 
be eccentric for ground-based (Advanced LIGO/ VIRGO 
and Einstein Telescope) and space-based detectors (DE- 
CIGO) 0. 

Including eccentricity in the waveform may be essen- 
tial for the detection of inspiraling eccentric binaries with 
matched filtering and to avoid a systematic bias in the 
parameter estimation |4lj . Using circular templates to 
detect waveforms with eccentricities eo > 0.1, leads to a 
significant loss of signal-to-noise ratio for grou nd-based 
detectors such as LIGO and VIRGO [II El]. A sim- 
ilar conclusion was reached for eccentric massive black 
hole binaries detected with LISA [Hj]. The orbital evo- 
lution and waveforms have been developed to first and 
second post-Newtonian (PN) order, including s pin -orbit 
and spin-spin contributions for eccentric orbits |45H50| . 

To assess the astrophysical impact of planned GW in- 
struments, it is essential to estimate the expected pa- 
rameter measurement precision of typical GW sources. 
This may be done by injecting a simulated GW signal 
into synthetic detector noise and carrying out a Monte 
Carlo Markov Chain (MCMC)-based matched filtering 
search for a parametrized template model to recover the 
posterior distribution function (PDF) of the estimated 
source parameters [Hlj]. Porter and Sesana [Hj investi- 
gated the cases of low-mass (100M Q ) and high (10 4 M Q ) 
mass black hole binaries on eccentric orbits using non- 
spinning, restricted 2 PN waveforms. They concluded 
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that eccentricity can significantly bias the recovered pa- 
rameters of the source for LISA if circular templates are 
used even if the eccentricit y is as small ase~ 10 -4 . More 
recently, Key and Cornish [52| extended that study by us- 
ing an effective 1.5PN waveform for inspiraling eccentric 
SMBHs [with m ~ (10 5 - 1O 7 )M ] taking into account 
eccentricity and spin effects in the template model. They 
found that the eccentricity measurement errors are of or- 
der Ae ~ 10 -3 for a range of mass ratios and a particular 
choice of angular parameters. 

Since the parameter space is large, 17-dimensional for 
an eccentric spinning binary, state-of-the-art MCMC cal- 
culations are numerically too expensive to explore the full 
range of source parameters. However, for a large signal- 
to- noise ratio (SNR), the PDF may be well approximated 
by an ellipsoid, and the parameter measurement errors 
can be estimated very efficiently using the Fisher matrix 
method [fl |4l[. Using this method, it has been shown 
that different source inclinations and sky locations lead 
to a wide range of parameter measurement errors sub- 
tending many orders of magnitude [9I4TTI Q [E3] . In this 
study, we carry out a Fisher matrix analysis to investi- 
gate the possible range of parameter estimation errors for 
eccentric binaries. 

Only a few studies have investigated the LISA param- 
eter estimation accuracy for eccentric inspiraling sources 
using the Fisher matrix method (cf. Refs. [p-fliL [53| 
for circular inspirals). Barack and Cutler [541 ] investi- 
gated the LISA errors for highly eccentric stellar mass 
compact objects inspiraling into a SMBH. They found 
that the influence of eccentricities on AAi/Ai ~ 1CP 4 
(error of the chirp mass), Aeo ~ 10~ 4 (error of initial 
eccentricity) and A^s ~ 10~ 4 (angular resolution error) 
is not substantial; the error estimates do not differ much 
from those obtained for circular orbits Q. However, they 
assumed only an arbitrarily chosen, single set of orienta- 
tions, which may not be representative of the typical er- 
rors. Yuncs ct al. [Hlf provided ready-to-use analytic ex- 
pressions for the Fourier waveforms of moderately eccen- 
tric sources. They have shown that eccentricity increases 
the detectable mass range of GW detectors toward higher 
masses by enhancing the orbital harmonics [HI, [l3[ • Yagi 
and Tanaka [HB] investigated the LISA errors for vari- 
ous alternative theories of gravity for spinning, small- 
eccentricity inspiraling SMBH binaries (eo ~ 0.01 at 1 
yr before merger), using restricted 2 PN waveforms, ne- 
glecting higher orbital harmonics and apsidal precession 
in the waveform. They have found that the eccentric- 
ity and the spin-orbit interaction reduce the parameter 
errors by an order of magnitude for spinning SMBHs in 
massive graviton theories, but not in Brans-Dickc-type 
theories. 

Neither of the previous systematic Fisher matrix stud- 
ies of parameter errors included the effects of relativis- 
tic pericenter precession for eccentric sources. However, 
precession effects introduce an additional feature in the 
waveform, and have the potential to break the degeneracy 
between parameter errors [l8j . In particular, spin-orbit 



precession has been shown to improve the source local- 
ization precision substantially during the last day of the 
inspiral [j| [ill EH ■ Similarly, GR pericenter precession 
may also be expected to improve the LISA parameter 
measurement accuracy. In fact, since pericenter preces- 
sion enters at a lower PN order, this improvement could 
take place well before the binary reaches merger. Lo- 
calizing the source before merger could be used to pro- 
vide triggers for electromagnetic (EM) facilities to search 
for the EM counterpart [19j . A coincident GW and EM 
observation of the same source could have far-reaching 
astrophysical implications [H, E3, Gil H3] 

In the present paper, we carry out a systematic pa- 
rameter estimation study for inspiraling SMBH bina- 
ries, taking into account both orbital eccentricity and 
the relativistic pericenter precession effect. We account 
for the evolution of the semimajor axis and eccentricity 
in our waveforms to leading order due to GW emission 
[42I [EH HH , but we neglect higher-order PN contribu- 
tions and spin effects. We compute the waveform in the 
frequency domain using the stationary phase approxima- 
tion (SPA, see Refs. [5H Io0l - [63l ]) and derive the signal- 
to-noise ratio (SNR) and the Fisher information matrix 
using a Fourier-Bcssel analysis for the parameter estima- 
tion of eccentric sources. To explore the possible range 
of parameter errors, we generate a Monte Carlo sample 
of binaries with random orientations and vary the masses 
and initial eccentricities systematically over a wide range 
relevant for LISA. We calculate the parameter errors for 
the standard three-arm LISA/NGO configuration, as well 
as for a descoped detector configuration, where one of the 
two independent interferometers is removed. 

In Sec. II, we summarize the basic formulas describing 
eccentric waveforms in the leading quadrupolc approxi- 
mation, using a Fourier-Bessel decomposition. In Sec. Ill, 
we derive the frequency domain waveforms and the LISA 
detector response. After a brief introduction of param- 
eter estimation using the Fisher matrix method in Sec. 
IV, we present results for specific systems in Sec. V. We 
summarize our conclusions in Sec VI. Some details of the 
calculations are described in Appendixes A and B. 

We use geometrical units G = c = 1. 



II. TIME-DEPENDENT ECCENTRIC 
WAVEFORMS 

To leading order, the waveform emitted by a binary 
moving on a Kcplerian orbit can be computed by the 
quadrupole approximation. In this approach the observer 
(i.e. the interferometric detector) is assumed to be far 
from the source and higher-order contributions; e.g., the 
effects of the spins and higher multipole moments are 
neglected, but the orbit is corrected for the effect of peri- 
center precession. For such precessing Keplerian orbits, 
the eccentric waveforms are given in Ref. [581 ] - We have 
rewritten the leading-order quadrupole tensor and trans- 
formed to the transverse-traceless gauge, which gives 
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Here is the true anomaly, which describes the azimuthal 
angle from the pericenter along the orbit as shown in Fig. 
[TJ The value 7 is the azimuthal angle of the pericenter 
relative to the coordinate system x axis in the orbital 
plane, e is the orbital eccentricity, a is the semimajor 
axis, Dl is the luminosity distance, is the inclination 
(the angle between the orbital plane and the line of sight 
to the observer), and m = mi+iri2, fJ- = mim^/m are the 
total and reduced masses (Figfl}. Using the well-known 
Fourier-Bessel decomposition, the polarization states can 
be expressed as a sum of harmonics of the orbital fre- 
quency [6l| 
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Here ft. = 4fim(aD is the amplitude, and = 
(S„ ± C n ) /2 and A n arc linear combinations of the 
Bessel functions of the first kind [J n (ne)] and their 
derivatives, 
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where a prime denotes the derivative, i.e. J' n {ne) 



n [J n _i(ne) + J n+ i(ne)] /2. 
Eqs. ([3}fl} are 



The phase functions in 
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nZ ± 2 7 , 
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(7) 



where I is the mean anomaly which is defined by the 
Kepler equation 



I = £ — esin£ = 2~nv(t — to) 



(8) 



In the Kepler equation £ is the eccentric anomaly, 
z/ = T _1 is the Keplcrian orbital frequency (here T = 
27tto -1 / 2 ci 3 / 2 is the Newtonian radial orbital period), and 
to is the time of pericenter passage (Hereafter, we set 



to = 0.) Equations (|6l7p show that the phase splits into 
a triplet due to the pericenter position 7. If the pericen- 
ter precesses, a triplet of frequencies appears in Fourier 
space for each harmonic [6l], [62[. Note that Eq. © is 
approximately valid during an orbit as long as w/c < 1 
and v = constant, but this equation requires modifica- 
tions on large time scales where the binary inspirals (see 
Eqs. [T2l - [T3l below), or at small separations, where the 1 
PN treatment breaks down. 

Pericenter precession leads to a time-dependent angle 
of pericenter, which may be written as 7(t) = 70 + 7(t) 
where 70 is the initial angle of the pericenter (Fig. [I]). 
Henceforth, we adopt pericenter precession from the clas- 
sical relativistic motion and assume the adiabatic evolu- 
tion of the orbital parameters. These effects are averaged 
over one radial oscillation period, i.e. (7) = A7/T, where 
A7 = 67rm[et(l — e 2 )] -1 is the angle of precession for an 
eccentric orbit governed by the geodesic equation of the 
Schwarzschild geometry (see e.g. Ref. \64\). In the fol- 
lowing we shall drop () for the average quantities, so wc 
write 

. _ 3m 3 / 2 _ 3?ti 2 / 3 {2nvf /3 

7 ~ a 5 /2(l_ e 2) - (1-e 2 ) ' (9) 

The 2.5 PN leading-order adiabatic evolution of the or- 
bital parameters due to gravitational radiation averaged 
over one radial period are [2l| 
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5tt(1 - e 2 ) 7 / 



24 



96 



3047W 5 / 3 (27rt/) 8 / 3 / 121 2 
G ~ 15(1 -e 2 ) 5 / 2 6 [} + 304 6 



[11) 



where M. = fi 3 ^ 5 m 2 ^ 5 is the chirp mass. (We used Ke- 
pler's third law, i.e. v = (2n)~ m 1 / 2 a~ 3 / 2 .) 

For an inspiraling system, the phase functions are 
= 2™ Jt^ v(t')dt' and < ± = ± 2 7o ± 

2 J t ^ oo A /(t')dt' . Equations ([6]), ([7]), are generalized as 
(here the t index is suppressed in Q n ±) 



Mt) v 
$„ = 2nn I —dv, 

J-oo v 

$„± = $„ ± 2 7o ± 2 / \dv 

. o v 



(12) 
(13) 
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FIG. 1. The geometry of an eccentric orbit. The coordinate 
system [x, y, z) is defined by the initial orbit, where the x axis 
points in the direction of the pericenter and the z axis is par- 
allel to the orbital angular momentum vector. In the reduced 
Kepler problem the body with mass /i = m\milm is orbit- 
ing the central mass m = m\ + 7712; the separation vector is 
r = eto(l — eo)/(l + eo cos 0), where eo is the orbital eccentric- 
ity; ao — m 1,/3 '(2ttuo) 2 ^ 3 (here uo is the orbital frequency) is 
the semimajor axis; <j> is the true anomaly (the angle between 
pericenter and the separation vector); and 70 is the pericen- 
ter position. The Kepler equation determines the evolution 
of the time parameter: £ — eo sin£ = 2nvo{t — to), where £ is 
the eccentric anomaly [tan£/2 = — eo)/(l + eo) tan <f>/2]- 
The adiabatic evolution of the eccentric orbit is driven by the 
pericenter precession (1 PN effect) and the inspiral (2.5 PN 
effect) of the compact binary due to gravitational radiation. 



e,v. 



m 



A 









K e ,v 


X 






— 1 ► 

y 



where <J>„± are phase functions which arise due to peri- 
center precession. Note that here one must incorporate 
the evolution in the eccentricity by solving Eqs. (jlOllip . 
i.e. v = v(y) = v\v, e(v)], and similarly for 7 [see Eq. (|36[) 
belowl . 



III. FOURIER TRANSFORMATION OF THE 
ECCENTRIC INSPIRAL WAVEFORM 

The sensitivity of a GW detector is usually given in 
Fourier space. Thus, to estimate the detection signal-to- 
noise ratio and measurement accuracy, we construct the 
Fourier transform of the waveform as 



phase approximation (SPA) can be utilized [62j (Ap- 
pendix B). We account for the adiabatic time evolution 
during the inspiral in the Fourier-transformed waveform 
h(f) using Eqs. an d the SPA. In this approxi- 

mation the Fourier transformation of the waveform be- 
comes a discrete sum over the harmonics of orbital fre- 
quency, /„ = nv. When the pericenter precession is taken 
into account, each harmonic /„, is split into a triplet 
f =(/„, f n ±) and therefore the waveform consists of the 
sum over these triplets of Fourier frequencies: 



h x (f) 



ho 



cos 6 



+S+A_ e J (*"- +7r / 4 >] , 
h+(f) = -^[si^e^Ae^"-^ 

n 

+ (l + cos 2 6) (s+A-e^*"— w / 4) 



(15) 



-B-A + e 4 (*"+-^ 4 ) 



(16) 



where /„ = nv, f„± = nv±\; h = 4M 5/3 (2tw) 2/3 /D l 
is the amplitude corresponding to the orbital frequency; 
and = 2irft n - $„, *, 1± = 2wft n ± - $„± are phase 
functions (where t n , t n ± are the time parameters of the 
SPA; see Appendix B). We have introduced the notation 

A± = (Twi^/Tr)" 1 / 2 and A = (nz>)~ 1/2 . The phases * n 
and $ n -t depend on the corresponding Fourier frequen- 
cies /„, f n ±, respectively. 

We recall that for circular orbits (i.e. e — >• 0) the wave- 
forms in Eqs. (|15I16[) simplify as 




^5/6^-7/6 

96 ir 2 / 3 D L 
5-^5/6^-7/6 



cos 6e l 



n 2 / 3 D L 



1 +cos 2 9) e 1 



(17) 
(18) 



where / = 2v is the (circular) Fourier frequency, and 
*± = 2nft c - $ c ± tt/4 + (3/4) (8TrMf)- 5/3 is the well- 
known phase function. 



A. LISA detector response 

With its three arms, LISA represents a pair of two or- 
thogonal arm detectors, / and //, producing two linearly 
independent signals. The frequency domain waveforms 
are 



Hf) 



h(t)e 



2-rritf 



(It 



(14) 



where / is the Fourier frequency. These integrals can- 
not be evaluated analytically without further assump- 
tions. However, since the orbital parameters (a, e) evolve 
very slowly relative to the GW phase, the stationary 



F$>"h x (f) 



Fl' H h + (f ) 
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with the antenna beam pattern functions 

= cos 2<fesin27/>s + Ms sin 20 s cos 2?As ,(20) 
■ cos 2<f>s cos 2ips — us sin 20s sm2ips , (21) 
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where hs,l = cos #s.l with (9s, (ps) being spherical an- 
gles of the source in the detector-based coordinate sys- 
tem. The angle ips is the polarization angle that can be 
expressed by the position of the detector and the orbital 
plane Q . The other antenna beam pattern functions are 
F+ x = x ((f>s — tt/4). The quantities 9s, (ps, and 
ips are time dependent, because the LISA constellation 
moves around the Sun, and these explicit time evolutions 
are @| 

Ms = 4f -^cos^ , (22) 
<Ps = ai (t) + £ + arctan ^{^f , (23) 
yj s = arctan cos^-cose(p s -V3A s cosfe) ^ 

where X s ,l = sin0 s ,L, PS,l = cos 6»s,l and </^ L = 
(p(t) — (ps.L, with 9s, (f>s being the spherical angles 
of the source's position. The angles 9l, <Pl corre- 
spond to the direction of orbital angular momentum 
in the bary center frame Jl] . In Eqs. (|22ll24j) . 6 = 
arccos [JIlPs + XlXs cos(4>l — (ps)] is the inclination [in 
Eqs. (|3l4p ] and the explicit time dependences are a\(t) — 
2nt/T - it/12 + a , (p(t) = (p + 2nt/T, and 

K = his. sin(0 L - fa) 

-^cos0(t) (fL L X s sin (ps ~ Ps^l sin<p L ) 

sin <p(t) (p, s X L cos (p L - p, L X s cos <p s ) .(25) 

We note that 9~l, <P~l are generally not constants for spin- 
ning binaries due to spin-orbit effects, but we neglect 
these effects here. 

We carry out the analysis for the single-detector case 
(J only) and the full two-detector configuration (I + II) . 

In practice, the measured signal in Eq. (|19[) is trun- 
cated at minimum and maximum frequencies correspond- 
ing to the start of the observation and the last stable orbit 
for each harmonic, respectively (see Sec. V below). 



IV. PARAMETER ESTIMATION 



In Eqs. (|27l28p p(n) is the probability for the noise, the 
inner product is defined by 

(9\k)=^f g -^ il df, (29) 

A;* is denotes complex conjugation and S(f) is the one- 
sided spectral noise density. The definition of the signal- 
to-noise ratio (SNR) of h is 

oo 

P 2 = (h\h) = mj h{f ^ f) df . (30) 



The waveform h(f) depends on the parameters A a which 
characterize the source. For a large SNR, the errors AA a 
have the Gaussian probability distribution 

p(AA c )= Po e- r " 6AAaAAV2 . (31) 

where po is the normalization factor and T a b is the Fisher 
information matrix defined by 

oo 

r ab = (d a h | db h) = mj ^mpil df , (32) 



with d a = d/dX a . The inverse of the Fisher matrix 
is approximately the £ a t, variance-covariance matrix for 
p> 1. which gives the accuracy of each parameter and is 
defined by S q6 = (r a6 ) _1 = (AA Q AA b ). The root-mean- 
square errors of the parameters A a are AA a = -\/S aa . 
For example, the error of the sky position solid angle 

is 

Afl s = 2n^ (Aji s Acp s ) 2 - (Aji s A^ s ) 2 . (33) 

The source localization sky area is an ellipse with 
semiminor and major axes (as,bs) given by Eq. (4.12) 
in Ref. [10]. The SNR and Fisher matrix for the LISA 
configuration are 

2 2,2 

P = Pi + Pn > 

r ab = T{ b + ril . (34) 

where the /, II subscripts distinguish the h 1 , h 11 wave- 
forms in Eq. ([T^|) . 



In this section we review the basics of Bayesian param- 
eter estimation. The measured signal s(t) is made up of 
the GW h(t) and the noise n(t) 

s(t) = h(t) + ri(t) . (26) 

We assume that the noise is stationary, Gaussian, and 
statistically independent at different frequencies. Then 
each Fourier component has a Gaussian probability dis- 
tribution and the different Fourier components of the 
noise are " uncorrected." i.e. 

p(n = n ) cx e - (nol " o)2 , (27) 
(n(f)n*(f')) = ±6(f-f')S(f). (28) 



V. MEASURING ECCENTRIC INSPIRALING 
SMBH BINARIES 

We focus on comparable-mass SMBH binaries in the 
range (10 4 — 10 7 )M Q , which corresponds to the measured 
frequency range 10 -4 to 10 _1 Hz. For initial configura- 
tions 1 yr before merger, we assume that the binary has 
orbital eccentricity eo and pericenter position 70. The 
ten-dimensional parameter space is 

A a = {lnD L , \nM,t c ,$ c , cp s , /2 S , <pL, PL, e , 70} 

In the circular case eo and 70 do not appear. Note that 
only one mass parameter, the chirp mass M., enters the 
leading-order waveform. Our assumptions are as follows: 
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To examine the effects of eccentricity and pcriccntcr 
precession, we neglect higher-order post-Newtonian 
(beyond 1 PN) orders and spins; we only use the 
heuristic pericenter precession in phase described 
above. 



In all cases, we take t c = $ c 
«0: 4>o = choice, as in Rcf. 



7o = 0. (We use the 
I) 



We assume that the observation time is 1 yr before 
the merger-more precisely, before the Newtonian 



last stable orbit (LSO), which is defined by [54 



; LSO 



1 



1 



'LSO 



27TTO I 6 + 2eLso 



3/2 



(35) 



where e^so is the final eccentricity at the last sta- 
ble orbit (p[e LS o) = vlso}- 

For the n orbital harmonic, the limits of inte- 
gration are taken to be f max = ^lso an d v m in = 
max{i/(i,/ c /n}, where Vq is the frequency 1 yr be- 
fore the LSO and f c = 0.03 mHz is the cutoff fre- 
quency of the LISA detector. 

We assume that the luminosity distance to the 
source is Dl = 6.4 Gpc, corresponding to a cos- 
mological rcdshift z = 1, and we use the comoving 
masses as free parameters, mf = (1 + z)rrii We 
do not take into account the Doppler phase due to 
the varying light travel during the LISA constella- 
tion's: orbit around the Sun. 

We parametrize the evolution of the orbital fre- 
quency with the instantaneous eccentricity follow- 
ing Ref. HH (Appendix A): 



v(e) = v 



*(e) 
c(e ) 



(36) 



where i>q and eo are the initial orbital frequency 
and eccentricity, and cr(e) follows from Ref. [21| . 



We truncate the harmonics at 
the signal power corresponds to [35 



> (l + eo) 
'(1-eo) 



1/2 
372 



where 99% of 



(37) 



where the bracket [J denotes the floor function (in- 
teger part of nonnegative argument). Here n max = 
{9, 24} for e = {0.3, 0.6}, respectively. 

We analyze 10 4 SMBH binaries where the angular 
variables were chosen randomly, i.e. for (j>g, <Pl 
in the range (0, 27r) and for 9$, @l in the range 
(-7r/2,7r/2). 



The computation of SNR and the Fisher matrix with 
the above general definition [Eq. Q3] is numerically ex- 
pensive for a large set of binaries. We resort to the SPA 
waveform. The SNR and the Fisher information matrix 
consist of three terms for each orbital harmonic which 
correspond to (/„, f n ±), respectively: 



•2 



p 2 = J2 (p 

n 

( f n 



Pn 



E 



' ah 



Pi-) 



■pn— 
1 ab 



(38) 
(39) 



where we have introduced the notations p 2 n n , n _ = 

(hn,n-\-,n— I ^n,n+,n— ); ^ab {,^a^n,n+ ,n— 

dbh n , n +,n-), and h n ,n+,n- = h(f n ,n+,n-)- Here we ne- 
glect the cross terms between different harmonics n, n+, 
and n— , in p and T a b- We use the LISA sensitivity curve 
generator [65J . In the SPA, we can change the integration 
variables from /„ ,f n ± to e: 



45R 



h n {e)h* n (e) ndv 



de , 



S [nv(e)] de 
d a h n {e)d b h* n (e) ndv 



S [nv(e)\ 



de 



de , 



(40) 



(41) 



where dvjde and v{e) are given by Eqs. (|A1IA3[) and 
e max = e LSO , e min = min{e c (n), e }. (Here e c (n) corre- 
sponds to f c /n, where f c = 0.03mHz is the cutoff fre- 
quency for the LISA detector.) 
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FIG. 2. (color online) Smooth probability density function of 
SNR for various initial eccentricities eo = 0.15,0.3,0.45,0.6 
and masses (10 6 — 10 6 ) Mq. The eccentricity dependence of 
SNR is almost negligible. 



VI. RESULTS AND DISCUSSION 

We have found that the LISA parameter estimation ac- 
curacy depends sensitively on the initial eccentricity and 
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TABLE I. The initial and final frequencies [vq and v\ — vlso) for various initial eccentricities (eo) and comoving masses 
(mi-rri2 with redshift 2 = 1) for a 1 yr inspiral before LSO. We use the shorthand notation ei = etso for the final eccentricity. 
We have completed with a dimensionless semimajor axis f — a/m at the initial (fo) and final points (fi). 
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FIG. 3. (color online). Smooth probability density function 
of SNR for various equal-mass binaries (for initial eccentric- 
ity e = 0.3). The SNR is 0(1O 2 ) for low-mass binaries 
(10 4 - 10 4 ) M Q . In the other cases, the SNR is C>(10 3 ). 



pericenter precession, and we have also examined the dis- 
tribution of parameter errors for a wide range of initial 
binary parameters and masses. The four angular param- 
eters (4>SiP'S, ( I ) l,P-l) are chosen randomly in a Monte 
Carlo sampling, and the cosmological redshift and lumi- 
nosity distance are fixed at z = 1 and Dj_, = 6.4Gpc. Fig- 
ures IMTUl show the histograms of the expected measure- 
ment errors of the binary parameters for the chirp mass 
AAi/Ai, initial eccentricity Aeo, and angular resolution 
Afls for equal-mass binaries with 10 6 M Q or 1O 7 M each. 
Our parametrization of the orbit is singular at eo = 0. To 
get around this, we use e = 10 -6 for circular orbits. We 
have presented three representative cases for the initial 
eccentricity: a nearly circular orbit with eo = 10 -6 (see 
Table HT1 and Fig. 0]), and orbits with medium eo = 0.3 



and high cq = 0.6 eccentricities. Our computations corre- 
spond to a 1 yr inspiral before LSO. The initial and final 
orbital frequencies (v$ and vlso) vary for the three kinds 
of initial eccentricities and different equal-mass SMBH 
binaries as shown in Table HI If the initial eccentricity 
eo increases, the initial frequency vq decreases 1 yr be- 
fore LSO, while the final frequency v^so does not change 
significantly, due to the fact that e^so is close to zero. 

Representative values are shown in Table |TT] for equal- 
mass SMBHs for a fixed set of angular configurations 
(05 = 4.642, fig = -0.3185, 4> L = 4.724 and p L = 
—0.3455). The table shows that accounting for the eccen- 
tricity in the waveform improves some of the parameter 
errors such as the errors of the angular resolution Afis, 
the initial eccentricity Aeo, and the chirp mass AAi/M 
for higher-mass SMBH binaries (10 6 - 10 7 )Af Q . For lower 
masses, i.e. (10 4 — 10 5 )M Q , the eccentricity and preces- 
sion have no essential effects on parameter estimation. 
For masses of 10 4 M Q , high eccentricity has no signif- 
icant effect on the parameters AM./M or Afig. How- 
ever, the initial eccentricity errors (Aeo) are improved for 
smaller masses typically by factors of 3-10 and they are 
greatly improved for larger initial eccentricities by orders 
of magnitude. Similarly, the source localization angu- 
lar resolution Afls decreases with increasing eccentricity 
and mass. However, pericenter precession does improve 
the parameter errors for higher-mass SMBHs. It can be 
seen that the eccentricity, compared to the circular or- 
bit case, does improve the error of luminosity distance 
ADl/Dl , but there is no essential change between the 
high and medium eccentricities with the inclusion of peri- 
center precession. The error of t c is not affected by the 
eccentricity or by pericenter precession. It is interesting 
to note that there are degeneracies (A<I> c ,A7o > 1) for 
errors of $ c and 70 in the nearly circular case, which 
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FIG. 4. (color online). Distribution of the major (top) and 
minor (bottom) axes (as, bs) of the sky position error el- 
lipse (Afis = nasbs) for various eccentric binaries with equal 
mass. (Here the pericenter precession is neglected.) The two 
panels correspond to 1 yr observation of (lO 7 — 10 7 ) Mq black 
hole binaries at z — 1 (Dl = 6.4Gpc) with LISA (2 detector). 
The angular resolution is improved for high-mass binaries. 



different parameter errors instead of an overall change in 
the SNR. 

Figure 0] shows the distribution of the major and/or 
minor axes of the sky position error ellipse for the nearly 
circular, medium and high initial eccentricity orbits. The 
shape of the error ellipse is important in coordinating 
GW observations with telescopes [TJ EH ■ It can be seen 
that the error of the ma jor and/or minor axes is improved 
for highly eccentric binaries. 
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can be explained by the fact that our parametrization 
of the orbit is singular at eo =0. For eccentric orbits 
(medium and high initial eccentricities) this degeneracy 
disappears(the errors of t c , $ c and 70 are not presented 
in Table HU). 

Figures [5] and [3] show the distribution of the SNR 
for different binary orientations, for various eccentricities 
and masses. The SNR is similar for equal-mass binaries 
with 1O 5 M < M < 10 7 M o , but significantly smaller 
for a SMBH of 10 4 Af Q or less. Remarkably, the SNR 
does not change significantly with the initial eccentricity, 
which is consistent with previous studies for small eccen- 
tricities [HI. This shows that the systematic improve- 
ment of the parameter estimation accuracy for eccentric 
sources is due to the breaking of correlations between 



FIG. 5. (color online). Estimated distribution of the chirp 
mass errors in the precessing and nonprecessing cases for 
the total (/ + 77, top) and single (I, bottom) detectors. 
The results are shown for medium (eo = 0.3) and high 
(eo = 0.6) initial eccentricities and higher-mass SMBH bina- 
ries (lO 6 — 10 6 ) Mq. For precessing sources the eo = 0.6 case 
is omitted in both figures due to the high degree of overlap 
with the eo = 0.3 case. 

Figures [5] and [JJ] show that the chirp mass errors 
are greatly improved for a larger initial eccentricity for 
10 6 M (;1 and 10 7 M Q equal -mass SMBH binaries (see also 
Ref. [131]). Furthermore, the chirp mass measurement 
errors are improved by an additional factor of 2-5 due 
to pericenter precession for relatively massive 1O 7 M0 bi- 
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FIG. 6. (color online). Same as Fig. [5] but for masses 
(10 7 - 10 7 ) M Q . 

naries, but not for 1O 6 M0 binaries. The typical chirp 
mass error is about 10 -5 for 10 7 M Q and 10 -4 for 1O 6 M 
binaries. 

Figures [7] and [8] show that the initial eccentricity er- 
rors are also improved for a high eccentricity, as the initial 
eccentricity parameter can be measured with high accu- 
racy; Aeo is about 10~ 5 to 10~ 4 for 10 7 M Q binariesm 
and about 10~ 4 to 10~ 3 for 10 6 M Q binaries. Pericenter 
precession improves the eccentricity errors by a factor of 
10 for 10 7 Af Q and by a factor of 2-3 for 10 6 M o . 

Figures l9l and ITOl show that the typical source sky local- 
ization accuracy for equal-mass binaries at z = 1 ranges 
between 10 -4 and 10 -2 steradians. Consistent with pre- 
vious studies [HI, EH , we find that the errors improve for 
higher initial eccentricities (eo = 0.6), compared to the 
cases of moderate to small initial eccentricities (eo = 0.3) 
for equal-mass 10 7 M Q binaries. The error Afig in the 
total two-detector case is about 1 order of magnitude 
better than for a single detector Q. For high initial ec- 
centricities, the angular resolution of the total detector 



FIG. 7. (color online). Estimated distribution of the ini- 
tial eccentricity errors in the precessing and nonprecessing 
cases for the total (I + II, top) and single (J, bottom) de- 
tectors. The results are shown for medium (eo = 0.3) and 
high (eo = 0.6) initial eccentricities and higher-mass SMBH 
binaries (l0 6 - 10 6 ) M . 



case is improved more compared to the single detector 
case for 10 7 Mq binaries (see Fig fTOj) . In contrast to the 
chirp mass and the eccentricity errors, the angular local- 
ization capabilities are not improved for eccentric equal- 
mass 10 6 M Q binaries but they are improved for 1O 7 M 
binaries. Figures IJJ] and [TU] clearly show that pericenter 
precession does not affect the sky position error for either 
mass choice. 

A possible explanation for the qualitatively different 
improvement of the sky position and mass-eccentricity 
errors is that the sky position is a slow parameter, as 
opposed to fast parameters like the chirp mass and ec- 
centricity [l8| . The slow parameters are determined by 
the slow orbital modulation of the signal by the detector's 
motion around the Sun, while the fast parameters also de- 
pend on the orbital phase. The correlations between the 
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FIG. 8. (color online). Same as Fig. [7] but for masses 
(10 7 - 10 7 ) M Q . 

slow parameters become large during the last week before 
merger when the signal-to- noise ratio increases, which 
prohibits the rapid improvement of the slow parameters' 
marginalized errors. Pcriccnter precession does not vary 
the binary inclination and cannot effectively break the 
correlation between slow parameters. However, percen- 
ter precession splits the GW frequency into a triplet for 
each harmonic which can break degeneracies for the fast 
parameters and efficiently improve their measurement er- 
rors. 



VII. CONCLUSIONS 

We carried out an extensive study of parameter esti- 
mation for eccentric binaries with arbitrary orbital ec- 
centricity. We computed the waveforms in the frequency 
domain by a new method optimized for taking into ac- 
count eccentricity, by changing the integration variable 
for the waveforms from the orbital frequency v(e) to the 



FIG. 9. (color online). Estimated distribution of the angular 
resolution AQ,s in the precessing case for the total (I + II, 
top) and single (/, bottom) detectors. The results are shown 
for medium (eo = 0.3) and high (eo = 0.6) initial eccentricities 
and higher-mass SMBH binaries (10 6 - 10 6 ) M . The curves 
for the nonprecessing eo = 0.3 and eo = 0.6 cases are omitted 
in both figures since they are close and very similar to the 
curves for the precessing ones. 



eccentricity variable e [38(. This results in an improve- 
ment of numerical precision as compared to standard ap- 
proaches in the frequency domain, where a Taylor series 
expansion of the orbital frequency v{e) (among others) 
in the eccentricity e is needed [55|. Our method is well 
suited for computing the Fisher matrix and the signal- 
to-noise ratio. Our parameter space is ten dimensional, 
consisting of four angles, the chirp mass, the luminosity 
distance, coalescence time and phase, initial eccentricity 
and pericenter position (compare Fig. [I}. The first eight 
parameters are standard for circular orbits too. 

We have examined the LISA parameter estimation er- 
rors for GWs emitted by eccentric inspiraling SMBH 
binaries including the effects of pericenter precession. 
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FIG. 10. (color online). Estimated distribution of the an- 
gular resolution AQs in the precessing and nonprecessing 
cases for the total (I + II, top) and single (/, bottom) detec- 
tors. The results are shown for medium (eo = 0.3) and high 
(eo = 0.6) initial eccentricities and higher-mass SMBH bina- 
ries (lO 7 — 10 7 ) Mq. For nonprecessing sources the eo = 0.6 
case is omitted in both figures since it is close and very similar 
to the curve for the eo = 0.6 precessing case. 



It is important to note that the angular resolution 
is significantly affected by the number of detectors (see 
Figs. 191 and [TTJ|) . However, nearly the same parameter 
estimation accuracy can be obtained for the single and 
total detector configurations for (10 6 — 10 6 )M Q binaries 
for fast parameters [181 ] like the chirp mass and eccentric- 
ity (Figs. [S]and[7]). The second detector systematically 
reduces the errors of these parameters for higher masses 
(10 7 - 1O 7 )M . 
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Appendix A: Orbital evolution and waveform 
According to Eqs. (|10lll[) . the equation 
18v 1 



du 
de 

can be solved as 



73 P 2 , 37 4 
24° ^96 



19 e(l 
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1 + 304 e 0, 
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v 1 - e o) ' 1S 

the integration constant that has been chosen to set the 
initial condition v(eo) = vq for the initial values eo and 
V Q . Then Eq. (|X2) is 



Based on a large set of simulated binary waveforms, we 
found that there is about 1 order of magnitude improve- 
ment compared to circular waveforms in LISA's angular 
resolution for highly eccentric sources (e.g. eo = 0.6) for 
relatively high SMBH masses ~ 1O 7 M . There is how- 
ever, a much smaller effect for lower-mass binaries in the 
range (10 4 — 10 5 )M Q . This improves the prospects for 
identifying the electromagnetic counterparts [TJ Ell of 
relatively high-mass eccentric SMBH mergers with LISA. 
Similar conclusions have been reached in Refs. (l2l. [l3|. 
However, we found that pericenter precession does not 
further improve the sky localization accuracy of the 
source, although it may further improve the measure- 
ment errors of mass and eccentricity parameters. 
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where c(e) = e 18 / 19 (l — e" 

From Eqs. (jlOllip one can compute the evolution of 
the time and phase functions [t — t c = ^fry, $ — 

$ c = 2tt J Q e -j^Tyrfe] in terms of eccentricity as [sec Eqs. 
(HE3)] 
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TABLE II. Parameter estimation errors for equal-mass SMBH binaries. The initial eccentricities eo are 10 6 (nearly circular), 
0.3, and 0.6. , the luminosity distance is Dl = 6.4Gpc [z — 1); and the angular parameters are cj)L = 4.724, ^ll = —0.3455 , 
<f>s = 4.642, and fi s = -0.3185. 
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where the I t and 1$ integrals are 
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with the constants a = 29/19, /3 = -1181/2299, 8 = 
-121/304, a = 11/1 9 and (3 = 124/2299. The inte- 
grals in Eqs. (|A6|A7j) can be evaluated with the Appcll 
functions which generalize the hypergeometric functions 
MM 

hie) - —£-Fi {—>f>>2>— (A8) 

19e 30 / 19 _ fa+1 ~ 5 + 3 . ,\ ,. n , 
U( e )= 30 2-F1 f ,(3,^-;8e 2 ) . (A9) 

To compute the time (AT) and phase (A«f>) differences 



the binary spends between the initial and final eccentric- 
ities eo and e\ during its evolution, Eqs. (|A4IA5[) are 
used, 

^-«^(^r '*<■»> <A10) 

Figures [11] and [T^] show the evolution of time and phase 
for various initial eccentricities, a fixed 1 yr inspiraling 
time before the LSO and 10 6 M Q equal-mass binaries. 
It can be seen that the eccentricity changes significantly 
near the coalescence, and the accumulated number of or- 
bits is decreasing for high initial eccentricity. 

Appendix B: Stationary Phase Approximation 

Consider the waveform h{t) = _4(i) cos $(i) with 
A(t)/A(t) < *(t) and $(f) < $(t) 2 (see e.g. Ref. ||), 
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FIG. 11. The evolution of the eccentricity as a function of 
time (as "lifetime" for the fixed 1 yr inspiraling time). The 
eccentricity changes significantly near the coalescence. 




FIG. 12. The evolution of the eccentricity in terms of the 
phase function for the fixed 1 yr inspiraling time. 



with its Fourier transform as 



F[A{t) cos $(t)] = / A(t) 



i*(t) 



i$(t) 



• 2 ™ tf dt . 



(Bl) 

To evaluate the Fourier integral one can use the 
stationary phase approximation (SPA). For an arbi- 

oo 

trary function of the time, / A(t)e m ^dt ~ 



^(7lA/2 7 r/*(7le i (* (r)+sisn [* (r) ] 7r/4 ), where the saddle 
point T satisfies <F(T) = 0. In Eq. (|B"T|) . the e 1 *^ terms 



have no contributions to the saddle point T . Moreover, 
= 2irtf — $>(t), and the stationary phase condition 
[*(T) = 0] implies that / = $(T)/(2tt). This provides 
a relation between the Fourier and orbital frequencies. 
Carrying out this exercise for an eccentric waveform con- 
sisting of many widely separated GW harmonics, the cor- 
responding Fourier frequencies are, respectively, f n = nv 
and f n ± = ni>± ^/n for the terms due to pericenter pre- 
cession. For circular orbits, the only nonvanishing term 
has frequency / = 1v. Therefore, the Fourier transform 
of harmonic functions with SPA are 

T[A{t)^{t)\ - ^^jmm^ r)1+f to 

^H(t)oo8*(t)] = «J^e^( ^"%3) 

where * [f(Vj] = 2?r f(T)t [v{T)] - $ [v{T)\ is the phase 
function and t [z/(T)], $ [v{T)] are derived from radiation 
reaction by Eqs. (|A4IA5p . 

Following Ref. [62| . the phase functions for eccentric 
compact binaries are 



= 2nft - $ n , 
= 2tt ft - , 



(B4) 
(B5) 



where the functions $„, <I>„± arc defined by Eqs. (|12I13[) 
and the first time derivatives are expressed as 



= 2nf - 2-rmv , 
*„± = 27r/ - 2-nnv =F 27 



(B6) 
(B7) 



There are three saddle points (t n , t n ±) following from the 
stationary phase conditions ^f n (t n ) = and ty n ±(t n ±) = 
0. It follows that there are three Fourier frequencies for 
each harmonic of the orbital frequency (denoted by /„, 
fn±)- The second time derivatives of the ^ n and 1 & n ± 
phase functions are 



ty n = —2-Knu , 
<F„± = —2-Knv =F 27 



(B8) 
(B9) 



where 7 is the time derivative of 7 induced by gravita- 
tional radiation; see Eqs. (|10lll[) . Then the phase func- 
tions of the waveforms, Eqs. (|15ll6p . can be expressed in 
terms of the time corresponding to the stationary phase 
and the acceleration of the pericenter precession, formally 

*«(/n) = 27T/ n i n (/ n ) - $„(/„) , (BIO) 
*«±(/n±) = 27r/ n± i„ ± (/„±) - $„±(/„±) . (Bll) 
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